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Chapter  1 

Finite  Difference  Schemes  for 
the  Boltzmann  Equation 


1.1  Introduction 

In  the  intermediate  report  [1],  Centered  Space  (CS)  finite  difference  schemes 
were  developed  to  solve  the  Lattice  Boltzmann  (LB)  evolution  equations  for 
the  distribution  functions  fi{r,t)  of  an  one  component  fluid  system,  when 
the  number  of  non  -  vanishing  discrete  velocities  was  N  —  6  (the  so  -  called 
seven  bit  model)  or  N  =  8  (the  so  -  called  nine  bit  model): 

fi(r,t  +  6t)  ~  fi{r,t)  -  —  [/i(r, t)  -  f*q(r,t)}  - 

Tp 

Ste i  ■  Vr/j(r,i)  +  -j^F(r,t)  ■  [es  -  u(r,t)]  /'’(r,  i) 

«  =  0,  1,  ...N  (1.1) 

In  the  above  equation,  as  well  as  in  the  whole  report  to  follow,  there  is  no  im¬ 
plicit  summation  over  the  index  i.  We  will  use  further  the  implicit  summation 
rule  (i.e.,  summation  over  repeated  indices)  only  for  cartesian  components 
denoted  by  Greek  letters  a,  /?,  7  . . while  summation  over  other  indices 
(including  i )  will  be  always  explicited,  when  such  operations  are  present. 

The  validity  of  the  CS  finite  difference  schemes  was  successfully  tested  for 
two  kinds  of  2D  flows  between  parallel  plates:  Poisseuille  flow  and  Couette 
flow.  The  correct  velocity  profiles  in  the  stationary  regime  (the  parabolic 


1 


one,  for  Poisseuille  flow,  and  the  linear  one,  for  Couette  flow),  as  well  as  the 
right  viscosity  value  were  recovered  for  N  =  6,  as  well  as  for  N  =  8  [1]. 

Our  early  attempts  to  use  the  CS  scheme  to  simulate  the  diffusion  phe¬ 
nomenon  in  a  two  component  system  revealed  an  unstable  behavior.  When 
this  scheme  is  used  to  get  the  time  dependent  solutions  of  the  lattice  Boltz¬ 
mann  equations  (1.1),  the  probability  functions  fi(x,  t)  become  negative  very 
often,  forcing  the  simulation  code  to  be  stopped.  Although  centered  finite 
difference  schemes  for  spatial  derivatives  arise  in  a  natural  way  and  thus, 
are  easy  to  introduce,  these  schemes  are  known  to  be  unconditionally  unsta¬ 
ble  in  accordance  to  the  von  Neumann  stability  analysis  [2,  3,  4]  applied  to 
hyperbolic  equations  like  the  Lattice  Boltzmann  equations  (1.1). 

This  chapter  reports  our  attempts  to  find  appropriate  finite  difference 
schemes  for  the  Lattice  Boltzmann  equations.  These  schemes  are  introduced 
for  the  nine  bit  model  using  the  characteristics  curves  of  hyperbolic  equations 
in  one  space  dimension  [4]  and  we  show  that  all  these  schemes  are  equivalent 
to  the  former  Lattice  Gas  Lattice  Boltzmann  (LGLB)  scheme  [1]  when  the 
Courant  -  Friedrichs  -  Lewy  number  CFL  (to  be  defined  further)  equals  unity. 
All  these  schemes  are  based  on  Lagrange  interpolation  procedures.  Some 
considerations  relative  to  the  treatment  of  boundary  lattice  nodes  are  also 
introduced  here. 


1.2  Spatial  derivatives  on  a  square  lattice 


To  solve  Eqs  (1.1)  on  a  square  lattice,  we  need  a  procedure  to  compute  the 
following  term  at  the  point  r  =  (x,  y) 


e*  •  Vr/i(r,  t)  =  (e*)aao/i(r,t)  (1.2) 

=  Vi  t)  4"  Ut  ^) 

(implicit  summation  over  cartesian  components  expressed  by  Greek  indices 
is  always  understood). 

A  possibility  is  to  use  the  Centered  Space  (CS)  finite  difference  scheme 


dxfi{x,y,t)  ~ 


fi(x  +  Sx,y,t)  -  fj(x  —  Sx,y,t) 
2  Sx 


(1.3) 


where  6x  is  the  lattice  spacing  in  the  x  direction.  A  similar  expression  holds 
also  for  the  derivative  in  the  y  direction  {dy);  the  lattice  spacings  are  all 
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equal  ( Sx  =  8y)  in  the  case  of  the  square  lattice.  As  mentioned  before, 
this  centered  finite  difference  scheme  was  found  to  be  inadequate  for  the 
simulation  of  diffusion  phenomena  in  a  two  component  system. 

Another  possibility  is  to  use  the  First  Order  Upwind  (FU)  finite  differ¬ 
ence  scheme  [3].  This  scheme  takes  into  account  the  fact  that  information 
propagates  along  the  vector  e*  and  thus,  one  should  consider  also  the  sign  of 
its  components  when  writing  the  corresponding  numerical  scheme.  For  the 
x  component,  we  have 


dxfi(x,y,t )  ~  < 


fi{x,y,t)  -  fj(x  —  8x,y,t) 
Sx 

fi(x  +  Sx,y,t)  -  fj{x,y,t) 
Sx 


5  (^i)x  ^  0 

(1.4) 


A  Second  Order  Upwind  (SU)  finite  difference  scheme  may  be  also  intro¬ 
duced  [3] 


(  3fj(x,  y,  t)  +  4 fj{x  +  5x,y,t)  -  f{(x  +  2Sx,y,t) 

2  8x 


dxfi{x,y,t )  ~  { 


(&i)x  ^  0 


3 fj{x,y,t)  -  4 fj(x-Sx,y,t)  +  fc(x  -  2Sx,y,t) 
2  Sx 


(6j)x  <  0 


(1.5) 

Although  the  FU  scheme  was  found  to  be  stable,  it  has  the  disadvantage 
of  a  lattice  spacing  dependence  of  the  viscosity.  This  is  a  general  feature  of 
first  order  schemes  (see  Section  1.10).  In  the  case  of  second  order  schemes 
(like  CS  and  SU),  no  lattice  spacing  dependence  of  the  viscosity  is  observed, 
but  these  schemes  are  not  stable  when  dealing  with  sharp  interfaces  (see 
next  chapter,  where  the  use  of  these  schemes  to  simulate  a  diffusion  couple 
is  discussed).  Other  schemes,  which  are  mainly  based  on  the  characteristics 
curve  and  use  an  interpolation  procedure,  are  described  below. 
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1.3  Characteristics 

We  start  from  the  differential  form  of  the  LB  equations  (1.1) 
dtfi{T,t)  +  e*  •  Vr/i(r,  t)  =  qi{r,t) 

i  =  0,  1,  ...  N  (1.6) 

where  the  source  term  is 

®M)  =  -- [/i(r,()  -  r(r,i)]  + 

Tp 

•  [e»  -  u(M)]  fi9{r,t) 

i  =  0,  1,  ...  N  (1.7) 

The  homogeneous  part  of  Eqs.  (1.6) 

=  dtfi(r,t)  +  e,  •  Vr/f(r,«)  =  0  (1.8) 

has  the  solution 

/i(r,  t)  =  fi{r0,t0)  =  constant  (1.9) 

since  the  vectors  e,  are  all  constant.  Consequently,  for  each  i  =  0,  1,  ...  N, 
the  current  point  moves  across  a  straight  line  (the  characteristics  line)  which 
passes  through  the  node  r  of  the  lattice  at  the  moment  t  and  has  the  orien¬ 
tation  of  the  vector  e*.  The  formal  solution  of  the  inhomogeneous  evolution 
equations  (1.6)  may  be  written  as 

fi(r ,t  +  St)  =  /,( r  -  e,«,  t)  -  -  /  [/i(r(«V)  -  ft* M«V)]  dt 

Tp  Jt 

+  tV  f +*  ■  tes  -  rwo.o 

k,B±  JSt 

i  =  0,  1,  ...  N  (1.10) 

where  the  integrals  are  calculated  along  the  characteristics  lines.  At  the 
initial  moment  t'  =  t,  the  current  point  is  r(t')  =  r  -  e^t  and  thus,  the 
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Figure  1.1:  Velocity  set  in  the  nine  bit  (N  =  8)  LB  model. 

integrals  in  the  equation  above  may  be  approximated  using  the  rectangle 
formula  to  get 

fi{ r,  t  +  St)  ~  fi{ r  -  ei8t,  t)  -  —  [  /*( r  -  etSt,  t )  -  /®9(r  -  e{St,  t )  ] 

Tp 

+  T~7?r  F(r  -  e<<5t,  t)  •  [  e*  -  u(r  -  e^t,  t)  ]  /Hr  -  t) 

kb± 

i  =  0,l,...N  (1.11) 

This  equation  provides  the  key  for  the  characteristics  -  based  finite  difference 
schemes  to  be  developed  further.  Eqs.  (1.11)  are  the  expression  of  the  well 
known  Lagrange  representation  in  classical  fluid  mechanics,  where  the  evo¬ 
lution  of  a  particle  on  its  trajectory  is  described  by  the  total  time  derivative. 

When  using  the  nine  bit  LGLB  scheme  on  a  square  lattice,  the  character¬ 
istics  lines  always  pass  through  the  lattice  nodes  at  any  time  step.  These  lines 
are  associated  to  each  velocity  speed  e*  in  Figure  1.1.  For  further  discussion, 
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Figure  1.2:  Lattice  points  along  the  characteristics  line. 


we  will  refer  to  one  of  these  characteristics  lines  (which  corresponds  to  any 
of  the  eight  velocities  e*  in  Figure  1.1).  We  may  choose  a  one  dimensional 
reference  system  whose  positive  direction  is  orientated  along  the  speed  ej  and 
we  represent  the  current  points  situated  along  this  line  as  in  figure  1.2.  In 
this  figure,  the  current  point  at  the  moment  t  is  denoted  by  j,  the  next  point 
on  the  characteristics  line  (which  becomes  the  current  point  at  the  moment 
t  +  St)  is  denoted  j  +  1,  while  the  preceding  one  is  denoted  j  —  1.  We 
should  remark  the  fact  that  the  point  j  (which  is  the  current  point  on  the 
characterisctics  line  at  the  moment  t )  always  refers  to  a  lattice  node,  but  the 
points  j  -  l  and  j  +  1  are  superposed  to  lattice  nodes  only  when  using  the 
LGLB  scheme,  as  seen  in  figure  1.3.  This  figure  shows  the  current  point  on 
the  characteristics  line,  at  different  time  steps,  in  the  X  -  t  plane,  when  the 
LGLB  scheme  is  used. 

When  Sx  is  the  spacing  between  the  latice  nodess  on  the  reference  axis  X 
of  the  characteristics  line  (note  that  Sx  may  be  equal  to  l  or  to  l  y/2,  where 


Figure  1.3:  Current  points  (o)  on  the  characteristics  line,  in  the  LGLB  model. 
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Figure  1.4:  Current  points  (o)  on  the  characteristics  line,  in  the  Finite  Dif¬ 
ference  LB  model  (CFL  <  1). 


I  is  the  side  length  of  the  elementary  square  lattice  cell)  we  have 

Sx  =  Xj  -  Xj- 1  =  Xj+ 1  —  Xj  (1-12) 

In  the  LGLB  model,  where  particles  are  considered  to  propagate  between 
lattice  nodes  during  one  time  step,  we  have 

Sx  =  cSt  (1.13) 


where  c  is  the  magnitude  of  the  propagation  speed  c  (c  €  {e*},  i  = 
1,  2  ...N)  and  St  is  the  time  step. 

When  the  lattice  spacing  Sx  does  no  more  equal  the  product  cSt,  the 
current  nodes  j  —  1  and  j  +  1  on  the  characteristics  line  passing  through 
the  lattice  node  Xj  (i.e.,  node  j)  are  no  more  superposed  to  lattice  nodes. 
These  current  points  lie  now  between  the  lattice  node  and  the  neighbors 
Xj- 1  or  Xj+ 1,  respectively.  This  general  case  is  shown  in  figure  1.4  and  is 
characterized  by  the  Courant  -  Friedrichs  -  Lewy  number 


CFL 


cSt 

Sx 


(1.14) 


When  CFL  =  1,  the  former  LGLB  case  is  recovered.  In  accordance  to  the 
general  theory  of  hyperbolic  equations  [2,  3,  4],  a  necessary  (but  not  sufficient) 
condition  to  be  satisfied  by  the  finite  difference  scheme  to  be  stable  is 


CFL  <  1 


(1.15) 
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1.4  Interpolation  procedures  on  the  charac¬ 
teristics  line 


The  Lattice  Boltzmann  equations  (1.11)  may  be  solved  using  an  iterative 
procedure.  The  value  /*( r  —  t)  of  the  distribution  function  (as  well  as 
the  value  u(r  -  erft,  t)  of  the  local  velocity)  may  be  calculated  using  an 
interpolation  procedure  along  the  characteristics  line  which  passes  through 
the  lattice  node  r  at  time  t  +  St.  The  current  point  at  the  moment  t  is 
r  —  eiSt,  which  corresponds  to  the  point  x  in  figure  1.5.  If  we  know  the 
values  of  the  distribution  function  at  at  the  moment  t  in  two  lattice  nodes  X\ 
and  X2,  we  may  use  the  first  order  Lagrange  interpolation  formula  to  evaluate 
the  distribution  function  in  the  point  x 


fiix,t)  =  7~r~r  Mx i»*)  +  *  _  *?  fi(x2,t)  (1-16) 

The  second  order  Lagrange  interpolation  formula  may  be  used  when  three 
values  of  the  distribution  functions  are  known 


fi{x,t)  = 


(x  -  x2){x  -  x3) 
(xi  -  x2)(xi  -  x3) 


fi(xi,t)  + 


(x  -  xi)(x  -  x3) 
(x2  -  xi)(x2  -  x3) 


fi(x2,t) 


+ 


(x  —  Xi)(x  -  x2) 
{x3  -  ®i)(*3  -  X2) 


fi(x3,t ) 


(1.17) 


Figure  1.5:  Current  points  (o)  on  the  characteristics  line. 
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Following  finite  difference  schemes  may  be  introduced,  depending  on  the 
choice  of  the  interpolation  formula,  as  well  as  on  the  number  and  the  location 
of  the  interpolation  points: 

1.  Upwind  scheme:  the  first  order  interpolation  formula  (1.16)  is  used 
with  x\  =  Xj~ i  and  X2  =  Xj. 

2.  Lax  Friedrichs  scheme:  the  first  order  interpolation  formula  (1.16)  is 
used  with  X\  =  Xj-\  and  x2  =  Xj+i- 

3.  Lax  Wendroff  scheme:  the  second  order  interpolation  formula  (1.17)  is 
used  with  X\  =  Xj- 1,  x2  =  Xj  and  x3  =  Xj+\. 

These  finite  difference  schemes,  as  well  as  the  space  centered  scheme,  are 
discussed  below. 

1.5  Upwind  scheme 

The  upwind  scheme  is  recovered  when  using  the  first  order  interpolation 
formula  (1.16)  with 


x  =  Xj  —  cSt 

X\  —  Xj — i 

x2  =  Xj  (1-18) 

In  this  case,  the  interpolation  formula  (1.16)  gives 

fi(xj  -  c5t,t)  =  fi(xj,t)  -  CFL  [fi(xj,t)  -  fi{xj-Ut)]  (1.19) 

A  particular  attention  should  be  given  to  boundary  nodes.  When  the 
vector  e*  is  pointing  outwards  the  lattice  domain,  the  interpolation  formula 
(1.19)  may  be  used.  But  the  interpolation  formula  should  be  somewhat 
different  when  th  vector  el  points  into  the  lattice  domain:  if  Xj  is  a  node 
located  on  the  boundary,  there  is  no  corresponding  Xj-i  node  in  this  case. 
The  node  Xj+ 1  may  be  used  instead,  and  we  have  the  following  extrapolation 
formula  in  this  case: 

fi(xj  -  c5t,t )  =  fi{xj,t)  -  CFL  [fi(xj+i,t)  -  fi(xj,t )]  (1.20) 


9 


After  interpolating  all  distribution  functions  fo,  i  =  0,  1  . . .  N  in  the  current 
point  r  -  e{5t,  the  local  velocity  u(r  -  eiSt,t)  should  be  evaluated  and 
introduced  in  Eq.  (1.11)  which  allows  to  calculate  the  new  values  of  the 
distribution  functions  in  all  lattice  nodes. 

Since  the  values  of  the  velocity  components  are  prescribed  on  the  bound¬ 
aries,  an  interpolation  procedure  becomes  necessary  to  evaluate  also  the  ve¬ 
locity  components  at  the  current  point  r  —  erft  when  the  node  r  is  located  on 
the  boundary.  Interpolation  and  extrapolation  formulae  similar  to  Eqs.  (1.19) 
and  (1.20)  give  the  value  of  the  velocity  u (xj  -  c8t ,  t),  which  allows  the  cal¬ 
culation  of  the  equilibrium  distribution  function  f%q(xj  -  c5t,  t)  in  Eq.  (1.11). 

A  very  important  feature  of  the  upwind  finite  difference  scheme  (1.19)  is 
the  recovery  of  the  former  Lattice  Gas  like  Latice  Boltzmann  (LGLB)  model 
when  CFL  =  1.  The  upwind  scheme  says  that  the  information  (i.e.,  the  dis¬ 
tribution  function  /*)  available  at  the  moment  t  at  the  current  point  Xj  —  c8t 
is  processed  there  as  a  result  of  collisions  (through  the  relaxation  term  mainly 
characterized  by  the  physical  relaxation  time  rp,  as  well  as  the  action  of  some 
force  F ),  and  thereafter  arrives  in  node  Xj  at  the  moment  t  +  5t.  But,  when 
CFL  =  1,  the  initial  current  point  Xj  —  c5t  is  identical  to  the  lattice  node 
Xj- 1,  in  accordance  to  Eq.  (1.19).  This  feature  belongs  also  to  the  Lax  - 
Friedrichs  and  the  Lax  -  Wendroff  schemes  to  be  discussed  below.  In  this  re¬ 
spect,  all  these  characteristics  curve  derived  finite  difference  schemes  provide 
a  natural  generalisation  of  the  former  LGLB  model  and  allow  to  separate 
the  lattice  spacing  Sx  from  the  particles  free  path  c5t.  This  separation  pro¬ 
cess  provides  the  keystone  to  the  LB  simulation  of  multicomponent  particle 
systems. 


1.6  Lax  -  Friedrichs  scheme 


The  Lax  -  Friedrichs  scheme  [2]  is  recovered  when  using  the  first  order  inter¬ 
polation  formula  (1.16)  with  the  following  interpolation  points 

x  =  Xj  —  cSt 

X\  —  Xj — i 

x2  =  Xj+ 1  (1.21) 


The  interpolation  formula  (1.16)  gives,  at  the  moment  t 


fi(xj  -  c6t,t) 


fi{Xj+\,t)  +  fj(Xj-i,t) 
2 
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CPL  fi{xj- jbjO 


(1.22) 


At  the  boundary  nodes,  the  Lax  -  Friedrichs  scheme  cannot  be  used,  but  the 
upwind  scheme  may  be  accepted  instead. 


1.7  Space  centered  scheme 

If  we  approximate 

+  fj(Xj- 1,^  _  (1.23) 

2 

in  the  Lax  Friedrichs  scheme  (1.22),  we  get  the  usual  space  centered  scheme 

Si{xs  -  cSt,  t)  =  /,(*„()  -  CFL  /<(*J+1’i)  ~  (1.24) 

This  scheme  is  known  the  be  unconditionally  unstable  [2,  3].  We  really  get 
negative  values  of  the  distribution  functions  when  using  this  scheme,  after 
a  certain  number  of  time  steps,  and  this  situation  forces  the  LB  code  to  be 
stopped.  The  number  of  time  steps  before  the  space  centered  computer  code 
crashes  is  dependent  on  the  CFL  number,  being  larger  when  this  number 
becomes  smaller.  Using  a  smaller  Courant  -  Friedrichs  -  Lewy  number  for  a 
given  lattice  spacing  Sx  is  unconvenient  since  this  means  a  smaller  time  step 
St,  which  generates  a  huge  increase  of  the  number  of  necessary  time  steps  to 
be  performed  in  order  to  recover  the  system  evolution  during  a  certain  time 
interval. 


1.8  Lax  -  Wendroff  scheme 

The  second  order  Lagrange  interpolation  formula  (1.17)  and  the  following 
interpolation  points 


x  =  x j  —  cSt 

Xi  =  Xj-i 
X2  =  Xj 
X$  Xj- (-1 
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(1.25) 
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1.9  Interpolation  Supplemented  Schemes 

1.9.1  General  Description 

One  can  see  that  the  procedure  described  by  Eq.  (1.11)  implies  the  following 
computation  steps: 

1.  Lagrange  interpolation  to  find  the  distribution  functions,  as  well  as  the 
local  velocities  in  the  points  r  —  e;<5i;  these  points  are  not  lattice  nodes 
and  are  all  different  for  i  =  1,  2,  ...  N 

2.  collisions  processing  and  the  addition  of  the  effect  of  external  forces 

3.  propagation  from  the  point  r  —  erft  to  the  lattice  node  r 

For  each  lattice  node  r,  one  has  to  perform  N  x  (N  + 1)  interpolations  since 
there  are  N  directions  e*  and  N  +  1  distribution  functions  to  be  interpolated. 
The  function  /0  should  also  be  interpolated  since  its  value  in  the  point  r  — 
is  needed  to  calculate  the  value  of  the  local  speed  u(r  —  e^i,  £);  the  local 
speed  u(r  —  et5t,  t)  is  necessary  to  calculate  the  equilibrium  distribution 
function  /,•( r  —  eiSt,t). 

An  equivalent  procedure  was  already  proposed  in  the  LB  litterature  [5, 
6,  7,  8].  This  Interpolation  Supplemented  Lattice  Boltzmann  (ISLB)  scheme 
(as  named  by  their  authors)  has  the  same  steps  as  described  above,  but  their 
order  is  changed: 

1.  the  collisions  and  the  action  of  external  forces  are  processed  in  the 
lattice  nodes  r 

2.  the  propagation  step  transports  the  result  towards  the  points  r  4- 

3.  the  interpolation  step  recovers  the  new  distribution  functions  in  the 
lattice  node  r  at  time  t  4-  5t 

To  avoid  confusion  with  another  scheme  to  be  described  below,  we  will 
denote  this  original  ISLB  scheme  as  ISLB-CP  (ISLB  Collision  -  Propagation) 
since  the  collision  step  is  done  before  the  propagation  step.  Figure  1.6  re¬ 
calls  the  main  characteristics  of  this  ISLB-CP  scheme.  In  this  figure,  lattice 
nodes  are  denoted  by  0,  1,  . . .  8.  Let  ro,  ri,  ...  r8  the  position  vectors  of  these 
lattice  nodes.  After  processing  the  effect  of  collisions  (as  well  as  the  effect 
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Figure  1.6:  The  ISLB-CP  scheme. 

of  external  forces)  in  the  original  lattice  nodes,  particle  distribution  func¬ 
tions  are  propagated  along  their  characteristics  lines  during  the  time  step 
St  and  arrive  in  the  displaced  nodes  O',  1',  . . .  8',  whose  position  vectors  are 
t'0,  r'1;  . .  .rg  (figure  1.6  shows  the  case  when  the  displacement  is  made  along 
the  e5  direction).  Thus,  for  each  vector  e*,  i  =  1,  2,  . .  .8,  the  collision  pro¬ 
cedure,  followed  by  the  propagation  procedure,  transforms  the  initial  set  of 
distribution  functions  {/j(rn,  t)},  n  =  0,  1  ...  8  into  the  set  {fi(r'n,  t  -I-  St)}, 
n'  =  O',  1',  ...  8'.  A  2D  interpolation  procedure  (to  be  described  further) 
allows  to  recover  the  value  of  fi(r0,t  +  St)  [5,  6,  7]  and  thus,  the  automaton 
algorithm  may  be  applied  again  for  the  next  time  step,  and  so  on. 

It  is  possible  to  imagine  an  ISLB  scheme  where  propagation  is  done  first, 
followed  by  collisions  processing.  This  scheme  is  presented  in  figure  1.7.  If 
we  start  from  the  Lattice  Boltzmann  Equations  (1.11)  we  may  write 

x  =  r  —  QiSt 

fi{x,t)  =  fi{r,t)  -  Stei-  Vr/i(r,i)  (1.29) 

and  thus 

fi(r,t  +  St)  ~  /i(x,t)  -  —  [fi{x,t)  -  f*q{x,t)}  + 

Tp 
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Figure  1.7:  The  ISLB-PC  scheme. 

F(x,f)  •  [es  —  u(x,  i)  ]  r(x,  t) 

i  —  0,  1,  . . .  N  (1.30) 

In  this  scheme,  which  we  will  denote  as  ISLB-PC  (since  the  Propagation 
step  is  done  before  the  Collision  step),  an  interpolation  procedure  in  two 
dimensions  is  used  to  determine  the  value  of  /i(x,  t)  in  Eq.  (1.30).  To  find  the 
value  of  /j(x,  t)  when  one  has  the  values  {fi{rn,t)},  n  =  0,  1  ...  8,  we  may- 
use  the  same  interpolation  procedure  as  for  the  ISLB-CP  scheme,  which  gives 
/i(ro,  t  +  St)  from  the  set  {fi(r'n,t  +  5t)},  nl  =  O',  1',  ...  8'.  The  ISLB-PC 
approach  may  be  viewed  as  being  equivalent  to  the  early  Finite  Difference 
schemes  developed  in  [1],  but  replaces  the  calculation  of  local  derivatives 
through  an  interpolation  procedure,  while  maintaining  the  propagation  of 
information  along  the  characteristics  line. 

The  ISLB-CP  and  ISLB-PC  schemes  use  biquadratic  interpolation  to  cal¬ 
culate  the  values  of  the  distribution  functions  in  the  central  node  (denoted 
by  0  in  Figures  1.6  and  1.7).  This  is  not  the  unique  possibility.  In  fact,  the 
authors  who  introduced  the  original  ISLB  scheme  [5,  6,  7,  8]  reccommmended 
the  use  of  an  upwind  biquadratic  interpolation  procedure.  This  procedure  is 
explained  briefly  in  Figure  1.8,  where  the  position  of  the  point  x  is  shown 
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Figure  1.8:  The  UISLB-PC  scheme. 
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in  two  cases,  when  the  characteristics  lines  are  directed  along  the  velocity 
vectors  e5  and  ei,  respectively  (the  other  velocity  directions  may  be  consid¬ 
ered  in  a  similar  way).  In  these  Upwind  ISLB-PC  (UISLB-PC)  cases,  the 
distribution  functions  are  interpolated  in  node  x  and  thereafter  propagated 
on  the  corresponding  characteristics  lines  towards  the  lattice  nodes  5  (Figure 
1.8a)  and  1  (Figure  1.8b).  An  Upwind  ISLB-CP  (UISLB-CP)  scheme  may 
be  introduced  in  the  same  way  as  the  UISLB-PC  scheme  described  here. 

The  Iterpolation  Supplemented  Lattice  Boltzmann  schemes  were  originally 
developed  to  deal  with  non  -  uniform  grids.  To  our  best  knowledge,  in  the  LB 
litterature  there  is  no  attempt  to  use  ISLB  schemes  to  multiple  component 
systems  whose  particles  carry  different  masses.  When  developing  the  ISLB- 
PC  scheme,  their  authors  started  from  the  early  LGLB  philosophy,  which 
considers  that  distribution  functions  are  redistributed  at  each  lattice  node 
because  of  collisions  and  are  moved  thereafter  in  the  direction  of  their  cor¬ 
responding  velocity  vectors.  When  compared  to  this  philosophy,  the  present 
discussion  of  Finite  Difference  Schemes  for  the  Boltzmann  Equation,  which 
is  based  on  the  characteristics  lines  of  the  Boltzmann  equation,  may  be  an¬ 
other  step  towards  a  more  rigorous  approach  to  Lattice  Boltzmann  models. 
Because  the  Lattice  Boltzmann  equation  is  a  hyperbolic  equation,  the  char¬ 
acteristics  curves  (which  are  straight  lines  since  the  velocity  vectors  et  are 
constant  [4])  may  provide  the  background  for  error  analysis  or  stability  in¬ 
vestigations. 

1.9.2  Interpolation  procedures  on  the  2D  lattice 

Both  ISLB-CP  and  ISLB-PC  schemes  need  an  interpolation  procedure  to 
evaluate  the  value  of  the  distribution  function  fi(x,t)}  at  the  point  x  = 
(x,y)  =  r0  -  eiSt  in  figure  1.7,  when  the  corresponding  values  {fi{rn,t)}, 
n  =  0,  1,  ...  8  of  the  distribution  function  are  known.  As  pointed  in  [8], 
second  order  interpolation  is  necessary  in  order  to  avoid  spurious  dependence 
of  the  fluid  viscosity  on  the  lattice  size. 

If  we  introduce  a  cartesian  coordinate  system  centered  in  the  node  r0  = 
{x0,yo)  in  figure  1.7,  such  that  the  X  axis  points  along  the  vector  ei  and 
we  use  the  fact  that  we  deal  with  a  square  lattice,  we  may  define  the  non  - 
dimensional  (natural)  coordinates  [9] 
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V 


(1.31) 


y  -  yo 

Sx 


G  —1,  +1 


To  simplify  the  notation,  in  this  subsection  we  will  write  /(x)  instead 
of  fi(x,t),  and  f(rn),  n  =  0,  1,  ...  8  instead  of  fi(rn,t).  The  interpolation 
procedure  gives 

/W  =  (1.32) 

n=0 

where  the  weight  coefficients  Nn  are  given  in  table  1.1  as  the  product  of  two 
second  order  Lagrange  interpolation  coefficients  (one  for  each  axis),  which 
are  defined  as  follows  (the  argument  6  stands  for  f  or  77) 


Li(0)  = 


2  fa\  - 


m 


[9-0][0-l] 
[(_!)  _  0][(-l)-l] 

[6>  -  (-1)][^  -  1] 

[o-(-i)][o-i] 

[g-(-l)][g-l] 


7-2  (Q\  _  l  V  /  1  L 

+[  >  ~  [1-  (-1)]  [1-0] 


(1.33) 


The  biquadratic  interpolation  procedure  described  above  may  be  applied 
to  any  lattice  node  which  is  not  situated  on  the  boundary  (i.e.,  bulk  nodes). 
Bulk  nodes  r0  have  all  the  eight  neighbors  rn,  n  =  0,  1,  ...  8.  When  the  node 
ro  is  located  on  the  lattice  boundary,  some  of  these  neighbors  are  missing, 
e.g.,  as  shown  in  figure  1.9,  which  refers  to  four  kinds  of  boundaries  (walls): 
bottom  (£),  top  (T),  left  (L)  and  right  (R),  Such  situations  are  encountered 
when  one  deals  with  flow  between  horizontal  or  parallel  plates  (we  will  not 
discuss  here  the  case  when  the  fluid  system  is  confined  into  a  box  -  this  case 
needs  to  consider  also  the  presence  of  corner  nodes). 

For  wall  nodes,  Eq.  (1.32)  changes  to 


/(x)  =  X)  Nnf(*n) 
n  6  A4v 


(1.34) 


where  the  sum  is  computed  with  respect  to  all  indices  in  the  set  J\fw,  W  G 
{B,T,L,R}.  The  elements  of  the  sets  Mb  are  (see  figure  1.9): 


Mb  =  {0,1,2,3,5,6} 
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Figure  1.9:  Lattice  nodes  on  the  walls:  (B)  -  bottom  wall;  (T)  -  top  wall; 
(L)  -  left  wall;  (R)  -  right  wall. 


Table  1.1:  Interpolation  weights  for  the  bulk  lattice  nodes. 


No 

mown) 

Ni 

luo 

n2 

n3 

Li(0L2M 

L&OLUn) 

n5 

LKOLKn) 

No 

Li{0L2+(v) 

n7 

Li  (OLlfo) 

N8 

Hit)  Lib) 

Mt  s  {0,1, 3, 4, 7, 8} 
Ml  =  {0,1, 2, 4, 5, 8} 
Mr  =  {0,2, 3, 4, 6, 7} 


We  may  use  linear  Lagrange  interpolation  for  the  Y  directions  when  deal¬ 
ing  with  horizontal  walls,  or  for  the  X  direction  when  dealing  with  vertical 
walls.  To  this  purpose,  we  define  the  corresponding  coefficients 


£»-(«)  = 


e  -  (-i) 

o  -  (-1) 


e  -  o 

(-1)  -  o 


L\+(0)  = 


e  -  i 
o  -  1 


£+•(»)  =  (1-35) 
while  the  interpolation  weights  Nn  are  given  in  table  1.2,  for  each  boundary. 
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Table  1.2:  Interpolation  weights  for  the  boundary  lattice  nodes. 


bottom 

N0 

Ll(OLh+W 

iVi 

L\(0  Ll+(n) 

n2 

N3 

L-(0  Lo+{v) 

n5 

N6 

L2M)L?(n) 

top 

N0 

m)L\-(n) 

Ni 

Ns 

Ll(£)Ll-(ri) 

l\(0  Ll-(n) 

n7 

LU()Ll~(v) 

N8 

£8(f  )£!■(<?) 

left 

No 

LUO  L2M 

Nx 

L\U  0  Ll(n) 

n2 

L\+(OL%(ri) 

n4 

Ll+(OLl(n) 

Ns 

LXUOLlin) 

Ns 

L\+(OLi(v) 

right 

No 

Ll-(0  Ll(n) 

n2 

L\-(OL%(n) 

Ns 

LIU 0  L2M 

Ni 

LlUOLiiO 

Ns 

L'UOLUn) 

n7 

L'siOLUn) 
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We  should  underline  the  fact  that  the  ISLB  schemes  (ISLB-CP,  ISLB-PC, 
UISLB-CP  and  UISLB-PC),  which  use  second  order  (biquadratic)  2D  inter¬ 
polation  procedures,  are  different  from  the  previously  introduced  Upwind, 
Lax  -  Friedrichs,  Space  Centered  and  Lax  -  Wendroff  schemes,  which  use  ID 
interpolation.  We  want  to  point  that  other  2D  interpolation  schemes  may 
be  also  considered,  besides  the  ISLB  schemes  (i.e.,  biquadratic  interpolation 
procedure  adopted  here  for  bulk  nodes,  associated  with  a  linear  -  quadratic 
procedure  for  boundary  nodes).  For  example,  the  bilinear  Lagrange  inter¬ 
polation,  as  well  as  the  linear  or  the  quadratic  elements  in  the  “serendipity 
family”  [9]  may  provide  some  alternatives. 

1.10  Numerical  simulations 

1.10.1  Computer  code 

In  the  remaining  section  of  this  chapter,  we  are  reporting  a  few  significant 
simulation  results  which  were  obtained  using  the  explicit  Finite  Difference 
Lattice  Boltzmann  (FDLB)  schemes  introduced  above.  These  results  refer  to 
2D  Poisseuille  flow  and  were  mainly  done  to  see  the  behavior  and  to  estimate 
the  errors  and  performances  of  each  scheme  we  investigated. 

To  reduce  the  programming  effort,  we  restricted  ourserves  to  the  nine 
bit  model  and  used  a  square  lattice.  The  computer  code  is  given  in  Ap¬ 
pendix  A.  All  the  necessary  simulation  parameters  (number  of  lattice  nodes 
in  the  X  and  Y  direction,  number  of  cycles  to  be  performed,  numerical 
scheme  to  be  used,  boundaries,  initial  configuration,  particle  masses,  local 
number  densities,  wall  velocities  and  so  on)  are  introduced  in  the  input  data 
file  wet  9.  input.  The  input  parameter  key  .scheme  controls  the  numerical 
scheme  to  be  used  during  the  simulation,  as  shown  in  Table  1.3.  For  further 
details  concerning  the  meaning  of  input  parameters  in  the  file  wet9 .  input 
please  refer  to  the  comment  lines  in  the  source  code. 

The  other  chapters  of  our  present  report  will  be  dedicated  to  the  physics  of 
diffusion,  surface  tension  and  wetting  phenomena,  as  it  can  be  recovered  using 
the  computer  codes  based  on  the  explicit  numerical  schemes  we  developed. 
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1.10.2  Single  Component  Fluid 

We  made  different  computer  runs  to  study  the  influence  of  the  parameters  Sx, 
8t  and  rp  on  the  kinematic  viscosity  v.  Initially,  the  flow  of  a  single  component 
fluid  within  a  2 D  channel  (Poisseuile  flow)  was  simulated  using  the  Upwind, 
Lax  -  Friedrichs  and  Lax  -  Wendroff  schemes.  A  5  x  25  lattice  was  used,  with 
the  same  forcing  term.  The  value  of  the  viscosity  was  determined  after  fitting 
the  parabolic  velocity  profile  across  the  channel,  as  done  in  our  intermediate 
report  [1]. 

Figure  1.10  shows  the  influence  of  the  lattice  spacing  5x  on  the  kinematic 
viscosity  v,  when  the  other  parameters  are  constant  (rp  =  10~* 1 2 3 * * * * 8  s,  St  =  10-9 
s,  particle  mass  m  =  1  amu,  temperature  T  =  300  K,  force  F  =  10-22  N). 
When  using  the  first  order  schemes  (Upwind,  Lax  -  Friedrichs  or  FU),  the 
viscosity  is  always  found  to  increase  when  the  lattice  spacing  8x  increases.  In 
particular,  when  using  the  Upwind  or  the  FU  scheme,  a  linear  dependence  of 


Table  1.3:  Possible  values  of  the  input  parameter  key_scheme,  which  controls 
the  numerical  scheme  to  be  used  during  computer  simulations. 


key_scheme  Numerical  Scheme  Section 


0  Upwind  (U)  1.5 

1  Lax  -  Friedrichs  (LF)  1.6 

2  Space  Centered  (SC)  1.7 

3  Lax  -  Wendroff  (LW)  1.8 

7  ISLB  -  CP  1.9 

8  UISLB  -  CP  1.9 

9  ISLB  -  PC  1.9 

10  UISLB  -  PC  1.9 

11  Linear  Serendip  Elements  1.9 

12  Bilinear  Interpolation  1.9 

13  First  Order  Upwind  (FU)  1.2 

14  Centered  Space  (CS)  1.2 


15  Second  Order  Upwind  (SU)  1.2 
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Figure  1.10:  Influence  of  the  lattice  spacing  Sx  on  the  kinematic  viscosity  u, 
as  determined  by  fitting  the  parabolic  velocity  profile  in  a  2D  channel:  (a)  - 
Upwind  scheme;  (b)  -  Lax  -  Friedrichs  scheme;  (c)  -  Lax  -  Wendroff  scheme; 
(d)  -  FU  and  CS  schemes. 
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Figure  1.10:  (continued)  Influence  of  the  lattice  spacing  5x  on  the  kinematic 
viscosity  v,  as  determined  by  fitting  the  parabolic  velocity  profile  in  a  2D 
channel:  (a)  -  Upwind  scheme;  (b)  -  Lax  -  Friedrichs  scheme;  (c)  -  Lax  - 
Wendroff  scheme;  (d)  -  FU  and  CS  schemes. 
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the  viscosity  on  the  lattice  spacing  Sx  is  found  when  this  parameter  becomes 
large  enough.  The  kinematic  viscosity  u  is  found  to  be  independent  on  the 
lattice  spacing  only  when  using  second  order  schemes  (Lax  -  Wendroff  or  CS). 
Similar  results  were  reported  when  using  the  ISLB  scheme  [5,  6].  A  rigorous 
analysis  of  an  ID  ISLB  model  [8],  revealed  that  one  should  use  a  second  order 
interpolation  formula  in  order  to  avoid  a  spurious  (lattice  spacing  dependent) 
viscosity  term  in  the  Navier  Stokes  equation. 

Despite  the  recovery  of  the  correct  (lattice  spacing  independent)  value 
of  the  viscosity,  the  Lax  Wendroff  scheme  in  our  computer  code  is  found  to 
become  unstable  when  the  Courant  -  Friedrichs  -  Lewy  number  exceeds  a 
certain  value  ( CFL  ~  0.04).  This  is  why  the  values  v  =  u(S: r)  in  figure 
1.10c  are  available  only  for  5x  >  7  x  10~5  m.  This  behavior  may  be  as¬ 
sociated  to  the  necessity  to  keep  a  small  value  of  the  Courant  -  Friedrichs  - 
Lewy  number  CFL,  in  order  to  preserve  the  validity  of  the  first  order  series 
expansion  (1.29),  which  is  essential  to  all  interpolation  based  FDLB  scheme. 
A  similar  problem  was  observed  when  using  the  CS  scheme  (Figure  1.1  Od) 
but,  in  this  case,  the  lattice  spacing  Sx  should  be  larger  than  1.5  x  10-5 
m  (which  means  CFL  ~  0.2)  in  order  to  avoid  the  negative  values  of  the 
equilibrium  distribution  functions. 

Figure  1.11  shows  the  influence  of  the  time  step  St  on  the  kinematic 
viscosity,  when  Sx  =  10-4  m  and  the  other  parameters  were  unchanged. 
When  using  the  Upwind,  Lax  -  Friedrichs  or  Lax  -  Wendroff  schemes,  the 
viscosity  v  is  always  found  to  decrease  when  increasing  the  time  step.  Only 
in  the  case  of  the  Lax  -  Wendroff  scheme,  the  dependence  v  =  v(St)  is  a 
linear  one,  in  accordance  to  the  general  theory  of  the  LGLB  model  [1,  10,  11] 

2 Tp  —  St  2  —  St  ksT  f  ^ 

u  =  — - x  c  =  — ~ - (1.36) 

2  2  m 

The  viscosity  is  found  to  be  independent  on  the  time  step  when  using  the 
Finite  Difference  Schemes  (FU  and  CS).  This  is  in  accordance  to  the  viscosity 
formula  derived  previously  in  our  intermediate  report  [1] 

v  =  Ttxc2  =  ^-kBT  (1.37) 

TTt 

Figure  1.12  shows  the  influence  of  the  relaxation  time  tp,  when  Sx  =  10-4 
m  and  St  =  10-9  s.  Eq.  (1.36),  which  predicts  a  linear  dependence  of  the 
viscosity  on  the  relaxation  time  tp,  is  found  again  to  be  valid  for  the  Lax  - 
Wendroff  scheme,  while  Eq.  (1.37)  is  well  verified  in  the  case  of  the  Centered 


Figure  1.11:  Influence  of  the  time  step  St  on  the  kinematic  viscosity  v,  as 
determined  by  fitting  the  parabolic  velocity  profile  in  a  2 D  channel:  (a)  - 
Upwind  scheme;  (b)  -  Lax  -  Friedrichs  scheme;  (c)  -  Lax  -  Wendroff  scheme; 
(d)  -  FU  and  CS  schemes. 
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Figure  1.11:  (continued)  Influence  of  the  time  step  8t  on  the  kinematic  viscos¬ 
ity  u,  as  determined  by  fitting  the  parabolic  velocity  profile  in  a  2D  channel: 
(a)  -  Upwind  scheme;  (b)  -  Lax  -  Friedrichs  scheme;  (c)  -  Lax  -  Wendroff 
scheme;  (d)  -  FU  and  CS  schemes. 
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Figure  1.12:  Influence  of  the  relaxation  time  rp  on  the  kinematic  viscosity  u, 
as  determined  by  fitting  the  parabolic  velocity  profile  in  a  2D  channel:  (a)  - 
Upwind  scheme;  (b)  -  Lax  -  Friedrichs  scheme;  (c)  -  Lax  -  Wendroff  scheme. 
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(d) 


Figure  1.12:  (continued)  Influence  of  the  relaxation  time  rp  on  the  kinematic 
viscosity  u,  as  determined  by  fitting  the  parabolic  velocity  profile  in  a  2D 
channel:  (a)  -  Upwind  scheme;  (b)  -  Lax  -  Friedrichs  scheme;  (c)  -  Lax  - 
Wendroff  scheme;  (d)  -  FU  and  CS  schemes. 
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Figure  1.13:  Dependence  of  the  kinematic  viscosity  v  on  the  inverse  of  the 
mass  (m-1)  of  the  particles  in  the  fluid  system,  as  determined  by  fitting  the 
parabolic  velocity  profile  in  a  2D  channel  (FU  and  CS  schemes). 

Space  (CS)  scheme.  Even  if  the  Upwind  and  the  First  Order  Upwind  (FU) 
schemes  give  also  a  linear  dependence  v  =  u(tp),  the  values  of  the  viscosity 
are  not  within  the  range  predicted  by  Eqs.  (1.36)  and  (1.37),  respectively. 

Although  the  Lax  -  Wendroff  scheme  gives  the  correct  dependence  of 
the  viscosity  on  the  simulation  parameters  Sx,  St  and  rp,  we  found  it  to 
be  unstable  when  dealing  with  a  homogeneous  two  component  fluid  system. 
It  is  this  reason  why  we  turned  later  to  other  schemes,  which  proved  good 
stability  for  such  systems. 

To  check  the  FU  and  CS  schemes  further,  we  performed  several  runs  with 
Sx  =  0.0001  m,  St  =  10-9  s,  rp  =  10-8  s  and  different  values  of  the  mass  m 
of  the  particles.  The  results  are  shown  in  Figure  1.13.  A  linear  dependence  of 
the  kinematic  viscosity  on  the  inverse  of  the  mass  (m-1),  which  is  predicted 
by  Eq.  (1.37)  is  found  only  in  the  case  of  the  CS  scheme. 

Table  1.4  and  Figure  (1.14)  show  the  dependence  of  the  kinematic  vis¬ 
cosity  v  on  the  lattice  spacing  Sx  when  centered  ISLB  schemes  are  used. 
Two  values  of  the  time  step  St  were  adopted,  while  the  relaxation  time  was 
maintained  constant  (tp  =  1  x  10~8  s).  One  may  see  that  table  1.4  has 
a  few  missing  entries  in  the  column  St  =  2  x  10-9  s.  These  missing  en- 
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Table  1.4:  Influence  of  the  lattice  spacing  Sx  on  the  viscosity  (ISLB  schemes). 


Sx 

[s] 

u 

[  m2/s  ] 

V 

[  m2/s  ] 

St  =  lx  10~9  s 

St  =  2  x  10-9  s 

ISLB-CP 

0.3  x  10“5 

0.020865 

0.5  x  10~5 

0.022264 

0.6  x  10"5 

0.022575 

0.021173 

0.7  x  10"5 

0.022785 

0.021425 

1.0  x  10~5 

0.023132 

0.021829 

2.0  x  10"5 

0.023457 

0.022198 

3.0  x  10~5 

0.023532 

0.022284 

4.0  x  10"5 

0.023560 

0.022317 

5.0  x  10~5 

0.023573 

0.022332 

10.0  x  10"5 

0.023592 

0.022354 

ISLB-PC 

0.3  x  10~5 

0.020059 

0.5  x  10"5 

0.021926 

0.6  x  10"5 

0.022332 

0.020750 

0.7  x  10"5 

0.022603 

0.021105 

1.0  x  10~5 

0.023093 

0.021665 

2.0  x  10"5 

0.023433 

0.022155 

3.0  x  10"5 

0.023521 

0.022265 

4.0  x  10~5 

0.023554 

0.022306 

5.0  x  lO"5 

0.023577 

0.022326 

10.0  x  HT5 

0.023592 

0.022352 
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tries  correspond  to  CFL  >  1,  when  the  numerical  scheme  is  always  unstable 
and  distribution  functions  become  negative.  For  both  St  =  lx  10~9  s  and 
St  =  2  x  10-9  s,  the  kinematic  viscosity  increases  when  increasing  the  lattice 
spacing  Sx  (i.e.,  when  decreasing  the  Courant  -  Friedrichs  -  Lewy  number  CFL 
below  unity).  When  the  lattice  spacing  is  large  enough  to  achieve  CFL  <  0.1 
(approximatively),  the  kinematic  viscosity  becomes  very  close  to  the  theoret¬ 
ical  value  calculated  in  accordance  to  Eq.  (1.36).  This  feature  is  seen  clearly 
in  Figure  1.14,  where  the  curves  v  =  u(6x)  exhibit  a  saturation  behavior. 
Thus,  the  error  in  the  determination  of  the  viscosity  becomes  negligible  (less 
that  one  percent  in  our  case)  when  the  first  order  series  expansion  (1.29) 
becomes  accurate  enough. 

In  order  to  see  the  influence  of  the  time  step  St  on  the  viscosity,  we  per¬ 
formed  several  computer  runs  with  the  ISLB  schemes,  using  a  fixed  value 
of  the  lattice  spacing  Sx.  This  value  was  choosen  to  be  Sx  =  3.0  x  10~5 
m,  which  lies  within  the  interval  where  the  approximation  (1.29)  is  accept¬ 
able.  The  results,  which  are  reproduced  in  Table  1.5,  as  well  as  in  Figure 
(1.15,  show  a  linear  dependence  u  =  u(St),  as  expected  in  accordance  to 
Eq.  (evisco).  This  equation  predicts  also  a  linear  dependence  of  the  kine¬ 
matic  viscosity  v  with  respect  to  the  relaxation  time  tp,  which  is  seen  in 
Figure  1.16,  for  Sx  =  3.0  x  10-5  m  and  St  =  1.0  x  10-9  s 

The  numerical  experiments  described  in  this  subsection  revealed  that  only 
the  CS  scheme,  as  well  as  the  ISLB-CP  and  ISLB-PC  schemes  give  the  correct 
dependence  of  the  fluid  viscosity  on  the  simulation  parameters,  for  a  single 
component  fluid. 
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Table  1.5:  Influence  of  the  time  step  St  on  the  kinematic  viscosity  v  (ISLB 
schemes). 


St 

V 

s] 

[  m2/s 

ISLB-CP 

0.5  x  10"9 
1.0  x  10“9 

1.5  x  10-9 
2.0  x  10-9 

2.5  x  10~9 
3.0  x  lO'9 

0.024165 

0.235323 

0.022906 

0.022284 

0.021665 

0.021047 

ISLB-PC 

0.5  x  10"9 

0.024159 

1.0  x  10~9 

0.023521 

1.5  x  10"9 

0.022891 

2.0  x  10“9 

0.022265 

2.5  x  10~9 

0.021642 

3.0  x  10~9 

0.021022 

(b) 


Figure  1.16:  Influence  of  the  relaxation  time  tp  on  the  kinematic  viscosity  v, 
as  determined  by  fitting  the  parabolic  velocity  profile  in  a  2D  channel,  for  a 
fixed  of  the  time  step  St  (1  x  10-9  s):  (a)  -  ISLB-CP  scheme;  (b)  -  ISLB-PC 
scheme. 
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1.10.3  Homogeneous  Two  Component  Fluid 

We  consider  a  homogeneous  two  component  fluid  system  subjected  to  Pois- 
seuille  flow,  whose  composition  may  be  varied.  The  system’s  composition 
may  be  described  using  the  mass  composition 


U) 


Pi 

Pi  +  P2 


Pi 

P 


(1.38) 


where  p\  and  p2  are  the  mass  densities  of  the  two  components  and  p  =  pi  +  p2 
is  the  mass  density  of  the  whole  system.  Another  possibility  to  account  for 
the  composition  of  the  homogeneous  system  is  to  use  the  mole  fraction 


ni  _  rii 

ni  +  n2  n 


(1.39) 


where  ni  and  n2  are  the  number  of  moles  belonging  to  each  component  and 
n  =  ni  +  n2  is  the  total  number  of  moles  in  the  system. 

In  order  to  study  the  composition  dependence  of  the  system’s  viscosity,  we 
kept  a  constant  number  of  moles  n  =  1  in  the  system  during  our  simulations, 
while  the  number  of  moles  belonging  to  the  first  component  was  varied  from 
0  to  1,  with  the  increment  0.1.  The  system’s  viscosity  was  determined  after 
fitting  the  parabolic  velocity  profile  established  across  the  channel,  as  done 
for  the  single  component  fluid. 

The  dynamic  viscosity  r)  of  the  system  should  equal  the  sum  of  the  dy¬ 
namic  viscosities  771  and  rj2  of  the  two  components 


V  =  Vi  +  m  (1-40) 

Since  the  dynamic  viscosity  77  of  a  fluid  is  expressed  as  the  product  of  its 
mass  density  p  and  kinematic  viscosity  v 


77  =  pv  (1.41) 

we  have,  for  the  two  component  fluid  system 

pv  =  {pi  +  p2)v  =  piVi  +  p2v2  (1-42) 

This  gives  a  linear  dependence  of  the  kinematic  viscosity  v  of  the  whole 
system  on  the  mass  composition  u 

v  =  —vi  +  —v2  =  u(vi~v2)  +  v2  (1-43) 

P  P 
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Thus,  the  slope  of  the  straight  line  v  =  v(u)  is  determined  solely  by  the 
difference  between  the  kinematic  viscosities  of  each  pure  component.  This 
slope  is  positive  when  vx  >  v2  and  negative  otherwise. 

Because  the  kinematic  viscosity  v  has  a  qualitatively  different  dependence 
on  the  relaxation  time  rp  when  using  the  Upwind  or  the  Lax  -  Friedrichs 
schemes  (see  Figure  1.12),  also  the  composition  dependence  of  the  viscosity 
of  a  binary  system  is  found  to  have  a  different  behavior  when  using  these 
schemes.  This  is  seen  in  Figure  1.17,  where  we  show  the  results  obtained  for  a 
binary  system  whose  particles  have  the  same  mass,  but  the  relaxation  times  tp 
are  different  for  each  component.  Although  the  linear  interpolation  schemes 
(Upwind  and  Lax  -  Friedrichs  schemes)  do  not  give  the  correct  dependence 
of  the  kinematic  viscosity  (1.36),  the  mass  concentration  dependence  of  the 
system  viscosity  is  always  linear,  as  predicted  by  Eq.  (1.43)  and  shown  in 
Figure  1.18. 

Figures  1.19  and  1.20  show  some  results  obtained  with  the  ISLB-CP 
scheme.  Similar  graphs  were  obtained  using  the  ISLB-PC  scheme.  The 
correct  dependence  of  the  viscosity  of  the  whole  system  on  the  mass  com¬ 
position,  as  well  as  the  other  simulation  parameters  (masses  mi  and  m2  of 
the  particles  belonfing  to  each  component,  time  step  St  as  well  as  relaxation 
times  t i  and  r2  of  each  component),  is  always  recovered. 

We  should  point  here  that  the  former  LGLB  model  cannot  account  for 
results  similar  to  those  shown  in  figure  1.19  for  a  homogeneous  two  com¬ 
ponent  system.  In  the  LGLB  model,  all  particles  share  the  same  thermal 
velocity,  even  if  their  masses  are  different.  Since  the  viscosity  in  the  LGLB 
model  is  mass  independent,  the  use  of  this  model  to  simulate  a  two  compo¬ 
nent  fluid  system  will  always  give  a  composition  independent  viscosity.  This 
erroneous  behavior  of  the  LGLB  model  is  limiting  seriously  its  application 
to  multicomponent  fluid  systems. 

To  compare  the  CS  and  the  ISLB-CP  schemes,  we  used  the  lattice  spac¬ 
ing  8x  =  0.0001  m  and  the  time  step  St  =  10~9  s  to  study  the  composition 
dependence  of  the  viscosity  of  a  very  particular  two  component  system.  The 
two  components  of  this  system  had  the  masses  mi  —  1  amu,  m2  =  2  amu 
and  the  relaxation  times  rx  =  1  x  10“8  s,  r2  =  2  x  10~8  s,  respectively.  In 
accordance  to  the  viscosity  formula  (1.37),  the  viscosities  of  the  two  compo¬ 
nent  fluids  are  identical  when  we  use  the  CS  scheme,  since,  for  this  particular 
system,  we  have 

T\  r2 


(1.44) 


(a) 


(b) 


Figure  1.17:  Composition  dependence  of  the  kinematic  viscosity  v  of  a  ho¬ 
mogeneous  two  component  fluid:  (a)  -  Upwind  scheme;  (b)  -  Lax  -  Friedrichs 
scheme.  The  two  components  have  the  same  mass  (1  amu)  but  the  relaxation 
times  are  different  (t\  =  10-8  s  and  72  =  2  x  10~8  s). 
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(a) 


mass  composition  p1/(p1+p2)  [%] 


(b) 


Figure  1.18:  Composition  dependence  of  the  kinematic  viscosity  v  of  a  homo¬ 
geneous  two  component  fluid  (Upwind  scheme).  The  two  components  have 
different  masses  (expressed  in  amu)  and/or  relaxation  times  (expressed  in 
s).  Curves  were  been  obtained  for  mi  =  1,  m2  =  0.5,  rx  =  1  x  10~8, 

r2  =  0.5  x  10~8,  (upper  curve),  mi  =  1,  m2  =  0.5,  tx  —  lx  10~8, 

r2  =  1  x  10-8,  (second  curve),  mx  =  1,  m2  =  2,  rx  —  1  x  10“8, 

r2  =  1  x  10-8,  (third  curve),  mi  =  1,  m2  =  2,  tx  =  1  x  10-8, 

r2  =  2  x  10~8,  (bottom  curve). 
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(a) 


mass  composition  p1/(p1+p2)  [%] 


(b) 


Figure  1.19:  Composition  dependence  of  the  kinematic  viscosity  v  of  a  homo¬ 
geneous  two  component  fluid  (ISLB-CP  scheme).  The  two  components  have 
different  masses  (expressed  in  amu)  and  the  same  relaxation  times  (expressed 
in  s).  Curves  were  been  obtained  for  m\  =  1,  m2  =  0.5,  T\  =  1  x  10-8, 
T2  =  1  x  10-8,  (upper  curve)  and  m\  =  1,  m2  =  2,  T\  =  1  x  10-8, 
r2  =  1  x  10-8,  (lower  curve),  respectively. 


42 


Figure  1.20:  Composition  dependence  of  the  kinematic  viscosity  v  of  a  ho¬ 
mogeneous  two  component  fluid  (ISLB-CP  scheme).  Curves  (from  top  to 
bottom)  were  been  obtained  for  (masses  are  expressed  in  amu,  while  relax¬ 
ation  times  are  expressed  in  seconds):  mi  =  1,  m2  =  1,  T\  =  1  x  10~8, 
r2  =  2  x  10-8,  mi  =  1,  m2  =  0.5,  T\  =  1  x  10-8,  r2  =  1  x  10-8,  mi  =  1, 
m2  =  2,  rl  =  1  x  10“ 8,  r2  =  2  x  10-8,  mx  =  1,  m2  =  1,  rl  =  1  x  10-8, 
r2  =  0.5  x  10~8,  mi  =  1,  m2  =  2,  ri  =  1  x  10-8,  r2  =  2  x  10-8. 
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Figure  1.21:  Composition  dependence  of  the  kinematic  viscosity  v  of  a  homo¬ 
geneous  two  component  fluid  (ISLB-CP  and  CS  schemes).  Curves  were  been 
obtained  for  mi  —  l  amu,  m2  =  2  amu,  r\  =  1  x  10"8  s,  r2  =  2  x  10-8  s. 


When  we  use  the  ISLB-CP  scheme  (as  well  as  the  ISLB-PC  scheme),  the 
viscosities  of  the  two  fluid  components  are  no  more  equal  since 


2ri  —  St 

mi 


2  r2  —  St 
m2 


(1.45) 


Thus,  the  viscosity  of  the  particular  homogeneous  two  component  system 
mentioned  above  is  composition  independent  when  the  CS  scheme  is  used, 
while  it  does  not  have  this  property  when  the  ISLB  schemes  (ISPB-CP  or 
ISLB-PC)  are  used  instead.  Figure  1.21  shows  the  computer  results. 
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Chapter  2 

Diffusion  Couple 


2.1  Description  of  the  model 

We  consider  a  two  component  system  whose  particles  carry  different  masses 
denoted  ma,  a  =  1,2.  The  system  is  described  by  two  distribution  functions 
sets  {  f°,N{ r,  t)  },  a  =  1,  2;  i  =  0,  1  . . .  N,  defined  on  a  2D  regular  lattice. 
The  distribution  function  ff’N(r,t)  gives  the  probability  to  find  in  node  r 
of  the  lattice,  at  time  t,  a  particle  of  mass  ma  having  the  velocity  ef,Ar, 
%  =  0,  1  ...  N.  As  discussed  previously  [1],  the  underlying  lattice  may  be  a 
square  one,  while  the  velocity  space  may  be  restricted  to  one  of  the  two  2D 
discrete  velocity  sets  {  e°’N  }  widely  used  in  LB  simulations,  which  correspond 
to  N  =  6  or  N  =  8.  These  velocity  sets  define  the  so-called  seven  bit  LB 
model  {N  =  6) 

©o’6  =  0  (2.1) 

<’6  =  ( eif  i  ei2  )  —  cos  27r~'6  ^  ,  sin  Ci  6 

i  =  1, . .  .6 

or  the  nine  bit  model  ( N  =  8) 

eo’8  =  0  (2.2) 


e?’8  =  cos 


7 r(i  —  1)  .  7r(®  —  1)  1 

2  ’  sm  2 -  C£r>8  ’  1 


=  1, . .  .4 
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,<t,8 


=  y/2 


7T  TtU  —  5) 

cos  1 T  +  2 


.  .  7 r  7r(z  -  5) 

■ sm  I  T  + 


C<7,8 


i  =  5, ...8 

The  “thermal”  velocities  ca^,  N  G  {6,8}  are  defined  by 


C<t,N  = 


kBT 

fTlaXN 


where 


Xn  = 


{;■ 


N  =  6 
N  =  8 


(2.3) 


(2.4) 


In  the  absence  of  external  forces,  the  distribution  functions  are  supposed 
to  satisfy  the  Boltzmann  equation  with  the  Bhatnagar  -  Gross  -  Krook  col¬ 
lision  term  ( dt  =  d/dt,  da  =  d/dxa ,  a  =  1,  2) 


(7  =  1,2  and  i  =  0,  1,  . . .  N  (2.5) 

and  therefore  evolve  in  accordance  to  the  Finite  Difference  Lattice  Boltzmann 
(FDLB)  equations 


f?N(r,t  +  6t)  ~  /"•>,«)  -  Ste’'N 


Ta 


[yr'W)  -  /r",e5 


a  =  1,  2  and  *  =  0,  1,  ...  N  (2.6) 

where  ra,  a  =  1,  2  are  physical  relaxation  times  and  f[’N,eq(r,t)  are  the 
equilibrium  distribution  functions. 

The  local  number  densities  are  defined  by 


N 


N 


n, 


*(r,t)  =  E  =  E  /r'V'e,(M) 


(2.7) 


i= 0 


i=0 
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while  the  local  density  and  the  local  velocity  of  each  component  are 
Ar(r,t)  =  mana(r,  t) 


N 


Mr>*)  = 


na{r,t)  S 

The  mass  averaged  (barycentric)  velocity  is  defined  by 
miniUi  +  m2n2U2  E£=i  P<*  u<^ 


(2.8) 

(2.9) 


u  = 


mini  +  m2n2  E£=i  P<? 


(2.10) 


The  mass  flux  jff  of  the  component  a  [12,  13]  is  defined  with  respect  to  the 
barycentric  velocity 

j<r  =  rhoa  (u,  -  u  (2.11) 

As  usual  in  LB  models,  the  equilibrium  density  functions  are  expressed 
as  power  series  in  the  components  of  the  local  equilibrium  velocity,  ua,N,eg 


fi 


a,N,eq  _ 


wfna{v,t) 

ea’N  .  n<r,N,eq  .  ucr,N,eq\2 

1  +  — - o -  +  V  * 


(2.12) 


(u<T,V,e?)2 

2  Xn<%, 


Xn<%,n 


where  the  weight  factors  wf  are 


_i_ 

2 

_1_ 

12 


for  the  seven  bit  model,  and 


4. 
9 

«f  =  1  i 

1 
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2  XncI 

i  =  0 
i  =  1, 

i  =  0 
%  =  1, 
i  =  5, 


N 


N  J 


,  6 


(2.13) 


4 

.8 


(2.14) 


for  the  nine  bit  model. 

Since  the  elements  of  the  velocity  sets  {e*,JV}  satisfy 


=  0 
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(2.15) 


E 


i 

Ecr.N  a,N  a.N  <r,N  ttr,N 
£ ia  ^i/3  ^i-f  ® iS 


X(T  Ctr,N  ^otfi 

0 

Xcr  C<r,N  (fiapfiyS  ~F  4 (S^aS  “I"  $8/3$ 7a) 


following  sums  are  easily 

i=cr 

f?,N,eq 

i=0 

computed 

=  na 

(2.16) 

i=a 

c ioL  Ji 
i= 0 

(2.17) 

i=<j 

ST  p<r  p<r  f<r>N>eq 
2s  e ia  ei/3 1  J  i 
i= 0 

=  [  Xcr  £(?  ^ot/3  "F  UaUp  j 

(2.18) 

i=a 

V'  pa'N f>a'N 

/  j  cia  cij3  c i'l  J  i 

i= 0 

=  Xa  4,N  (  (5a/3  ^7  T  <5a7 Up  +  S/3yUa  ) 

(2.19) 

As  discussed  previously  [1],  the  macroscopic  behavior  of  a  single  com¬ 
ponent  fluid  system  is  independent  of  the  discretization  of  the  phase  space 
(i.e.,  independent  of  the  number  of  elements  N  of  the  discrete  velocity  set 
{ef,Ar}  used  in  the  FDLB  model).  In  the  case  of  a  two  component  system, 
we  will  consider  the  same  model  (seven  bit  ot  nine  bit,  i.e.,  the  same  N )  for 
both  components.  The  choice  N  =  6  or  N  =  8  is  only  a  technical  detail 
related  to  the  FDLB  computer  code  and  has  no  influence  on  the  results  at  the 
macroscopic  scale.  To  reduce  the  number  of  indices  in  the  following  formulae 
in  this  chapter  (as  well  as  in  the  subsequent  chapters  of  this  report),  we  will 
discard  the  index  N.  Consequently,  we  will  write 


dtff  (r,t)  +  4A/fM)  = 


a  —  1,  2  and 


instead  of  Eq.  (2.5), 
fi’eq  =  Wi  na{r,t) 


1  + 


e?  •  \ia,eq 
X*4 


+ 


—  E/f(r.O  - 

fi’eq(r,t)} 

To- 

*  =  0,  1,  ... 

N 

(2.20) 

(ef  •  u<r’e9)2 

WA 

(u<T,e9)2  - 

2  x*4 

(2.21) 

instead  of  Eq.  (2.13),  and  so  on. 
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2.2  Chapman  -  Enskog  expansion 

The  distribution  functions  are  usually  expanded  as  power  series  in  the  Knud- 
sen  number  e 

It  =  ftm  +  e/r(l)  +  eVf'2'  +  •  •  •  (2-22) 

while  two  corresponding  time  scales  and  one  length  scale  are  introduced  when 
computing  time  and  space  derivatives 

dt  =  edh  +  £20(1  (2.23) 

d„  =  eda  (2.24) 


After  introducing  the  Chapman  -  Enskog  expansion  (2.22),  as  well  as  the 
corresponding  expansions  (2.23)  and  (2.24)  of  the  time  and  space  derivatives 
in  the  Boltzmann  equation  (2.20) 

(edh  +  e2db)[ftm  +  ef?m  +  e2ftm  +  ■■■' 

l  + 

A[f’m  +  eftm  +e2ft(2)  +  ...' 

l  - 

-  —  [/”<0)  +  <=//(1)  +  e2ftm  +  ■■■-  ft’"} 

ra 

(2.25) 

we  can  separate  the  zero-th,  first  and  second  order  Boltzmann  equations  with 
respect  to  the  Knudsen  number  e 

o 

H  £ 

i 

ii 

o 

-  f[’eq]  (2.26) 

d,Jtm  +  deftme 

•i a 

(2.27) 

d,2ftm  +  dtJtW  + 

•a 

(2.28) 

From  the  zero-th  order  Boltzmann  equation  (2.26),  we  get,  for  any  o  e  { 1,2} 
and  i  =  0,  1  ...  N 

ftm  =  Stm  (2-29) 

This  condition  ensures  also  the  validity  of  the  zero-th  order  mass  and  momen¬ 
tum  equations  (2.32)  and  (2.35),  which  will  be  discussed  in  the  next  section. 
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Taking  into  account  the  expression  (2.7)  of  the  local  number  density  na  and 
the  series  expansion  (2.22),  we  get 

i=N 

£  f°{l)  =0  V/  >  1  (2.30) 

i=0 

As  usual  in  many  Lattice  Boltzmann  models  for  multicomponent  fluids  [13, 
14,  15,  16,  17],  we  assume  that  all  components  have  the  same  equilibrium 
velocity  in  the  absence  of  external  forces  and  long  range  interactions 

<’e<?  =  u'Q  V<r  =  1,2  and  a  =  1,2  (2.31) 

Consequently,  the  equilibrium  distribution  functions  f°'eq  =  /,•  are  ex¬ 
pressed  as  series  expansions  (2.21)  with  respect  to  the  components  u'a  of  this 
equilibrium  velocity.  At  this  stage,  we  are  not  interested  in  the  procedure 
which  should  provide  the  means  to  compute  the  values  of  the  components 
u'a ;  this  procedure  will  be  discussed  later. 

2.3  Conservation  equations 

The  zero-th,  first  and  second  order  mass  conservation  equations  for  each 
component  a  =  1,  2  are  recovered  from  Eqs.  (2.26  -  2.28)  after  summation 
with  respect  to  i  (multiplication  by  ma  is  omitted) 

-  —  E  [/f"”  -  =  0  (2.32) 

T°  i=0 

dh  E  f°{0)  +  dp  E  f’^elp  =  -  T  E  f!m  (2.33) 

i=0  i—0  Tcr  1=0 

/r(0)  +  Sh  E  +  9p  E 

2=0  *=0  i=0 

The  zero-th,  first  and  second  order  momentum  conservation  equations  are 
also  recovered  from  Eqs.  (2.26  -  2.28)  after  multiplication  with  and 

summation  with  respect  to  i  and  a 

-  e  —  e  [/r<0)  -  sr  ]  c  =  o  (2.35) 

cr=l  i=0 


1  i=N 

-f  E/” 

'<r  i=0 


(2.34) 
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<7=2 

dtl  53  m' 

i=N 

v  fcr(0V<T  4- 

<7  Z^l  Jl  lOi  — 

(7=2 

dp  53 m 

i=N 

V  fa{Q)ea  eaa  - 

o'  /  >  n  '-’icr'ip 

(7=1 

i—0 

<7=1 

i=0 

(7=2 

™  i=N 

-  £  - 

^  E  /r(1)<£ 

(2.36) 

<7=1 

i=0 

(7=2 

i=N 

<7=2 

i=N 

a,  E 

53  fi{  )eia 

+  53 

m*  53  /r(1)<a  + 

(7=1 

i=0 

(7=1 

i=0 

(7=2 

i=-N 

<7=2 

m  i=N 

dp  53  ma 

(7=1 

E 

t=0 

=  -£- 
<7=1 

/  v  J  i  c ia 

t=0 

(2.37) 

The  zero-th  order  mass  and  momentum  equations  (2.32)  and  (2.35)  are 
automatically  satisfied  in  accordance  to  Eq.  (2.29).  Moreover,  the  right  hand 
sides  of  the  first  and  second  order  mass  equations  (2.33)  and  (2.34)  vanish 
because  of  Eq.  (2.30)  and  thus,  these  equations  are  rewritten  as 

0tl  E  /f<0)  +  Sg  E  /f‘%  =  0  (2-38) 

i= 0  t=0 

dH  E  /r(0)  +  E  /f(1)  +  ^  E  =  0  (2.39) 

i=0  i=0  i=0 

2.4  Mass  equations 

Since  =  f[’eq,  we  can  use  the  series  expansion  (2.21),  as  well  as  the 
properties  (2.16)  and  (2.17),  to  rewrite  the  first  order  mass  conservation 
equation  (2.38)  in  a  more  familiar  form,  after  multiplication  by  ma 

dtlPa  +  9p(pau'p)  =  0  (2.40) 

Here 

Pa  =  manff  (2.41) 

is  the  local  density  of  the  component  a. 

Using  the  first  order  Boltzmann  equation  (2.27)  to  express  in  the 
second  order  equation  (2.28),  we  get 

a,J°m  +  -  27V0„  +  0„/f<%]  - 
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r^dff^ele^  =  (2.42) 

•O’ 

This  equation,  multiplied  by  ma  and  summed  over  i,  gives  the  following  form 
of  the  the  second  order  mass  conservation  equation  (2.39) 

dt2pa  +  ra{dh)2pa  ~  Tadadp (nakBT 5ap  +  Pau'au'p)  =  0  (2.43) 

To  derive  this  result,  we  used  the  first  order  mass  equation  (2.40),  as  well  as 
the  property  (2.18)  and  the  definition  (2.3). 

To  recover  the  mass  conservation  equation  for  the  component  a  up  to 
the  second  order  with  respect  to  the  Knudsen  number,  we  sum  together 
Eqs.  (2.40)  and  (2.43)  multiplied  by  s  and  e2,  respectively,  and  take  into 
account  the  expressions  (2.23)  and  (2.24)  of  the  time  and  space  derivatives 

dtPa  +  da(pau'a )  +  Ta{dt)2Pc  - 

ra  8apdadppa  -  Tadadp(pau'au'p)  =  0  (2.44) 

ma 

When  summing  over  a, 
system 


a= 2 

kBT  dadp 

<7=1 

where 

P  =  E  P „  (2.46) 

<7=1 

is  the  local  density  of  the  whole  fluid.  In  the  particular  case  ra  =  r  Va,  the 
mass  equation  for  the  whole  fluid  (2.45)  becomes 

dtp  +  da{pu'a)  +  r(dt)2 pff  -  t dadp p  -  rdadp  (J pu'au'p )  =  0  (2.47) 

where 

P  =  kBT  E  n,  (2.48) 

<7=1 


we  get  also  the  mass  equation  for  the  whole  fluid 


<7=2 


8tp  +  da{pu'a)  +  ( 8t )2  J]  Tcjpa  ~ 


a— 1 


Taper 

ma 


-  dadp 


^  T<rP"  ^  U<*UP 


=  0  (2.45) 
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is  the  local  pressure  of  the  fluid  system,  which  is  a  mixture  of  two  ideal  gases. 

The  mass  equation  for  a  single  component  (2.44),  as  well  as  the  mass 
equation  for  the  whole  system  (2.45)  contain  the  first  and  second  order  time 
derivatives  of  the  local  density.  The  presence  of  both  time  derivatives  is  a 
characteristics  of  the  telegraphist  equation  [18],  which  describes  a  propaga¬ 
tion  phenomenon.  Consequently,  one  may  expect  that  the  density  profile 
across  the  system  exhibits  a  kink  (wrinkle  or  peak),  i.e.,  a  non  monotonic 
propagating  pulse.  This  behavior  should  be  dominant  for  small  values  of  the 
product  between  the  local  density  and  the  equilibrium  velocity  {pau'a  ~  0), 
when  the  mass  equation  for  the  component  a  reduces  to  the  true  telegraphist 
equation 

9,/V  +  =  2>VV  (2.49) 

where  the  diffusion  coefficient  is 


mff 


(2.50) 


and 


(2.51) 


is  the  finite  speed  at  which  information  travels  in  the  system.  The  diffusion 
equation  is  recovered  in  the  case 


V 


(2.52) 


2.5  Equilibrium  velocity  u' 

The  right  hand  side  of  the  Boltzmann  equation  (2.20)  represents  the  collision 
term  which  expresses  the  variation  of  the  equilibrium  distribution  function 
as  a  result  of  collisions  between  particles.  Collisions  should  preserve  the  local 
momentum  of  the  whole  system.  Since  the  momentum  equation  is  derived 
from  the  Boltzmann  equation  after  multiplication  by  mffefQ  and  summation 
over  the  indices  a  and  i,  we  should  have 

[/f  -  x~]  4.  =  0  (2-53) 

<7=1  T<r  i=  1 
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Since  f°'eq  is  expressed  as  a  series  expansion  (2.21  in  the  equilibrium  velocity 
u9a  (which  is  supposed  to  be  the  same  for  all  a),  we  can  use  Eq.  (2.17)  to  get 


& — 2  njj  rp\  (J—-2  yu-  j—N 

/  ut'a'i'cr  ~<r 

««  E  — —  =  E  —  E  fi  eia 

<7=1  'ff  <7=1  t=l 


(2.54) 


which  allows  the  determination  of  the  components  u'a  of  the  equilibrium 
velocity  in  the  node  r  of  the  lattice. 

If  we  take  into  account  the  fact  that  the  distribution  functions  ff  are 
expressed  as  power  series  (2.22)  in  the  Knudsen  number,  as  well  as  the  fact 
that  =  fi’eq,  we  get  the  following  relations  from  Eq.  (2.54) 


771  *  ^  n\ 

E  —  E  if><&  = «  .«>i  (2-55) 

<7=1  T<*  i=l 


2.6  Momentum  equation 

If  we  use  Eq.  (2.55),  the  first  order  momentum  equation  (2.36)  becomes 


<7=2  i=N 


<7=2  i=N 


dt,  E  E  +  SI>  E  E  =  0  (2-56) 


<7=1  i=0 


<7=1  t= 0 


Introducing  the  relations  (2.17)  and  (2.18)  in  the  above  equation,  we  get  the 
Euler  equation 


<7=2  <7=2 

dtl  E  P°U'a  +  d<*P  +  dp  E  P<X«/3  =  0 

<7=1  <7=1 


(2.57) 


Using  the  first  order  Boltzmann  equation  (2.27)  to  express  in  the 
second  order  momentum  equation  (2.37),  we  get 


<7=2  i=N 


<7=2  i=N 


dt2  E  E  fi{0)eia  +  iptif  E  T°m°  E  /r(0)<a  - 


<7=1  i=0 


<7=2  <7=2 

2<9tl  dh  E  P<X  +  5aP  +  ^E  /V«au0  - 

<7=1  <7=1 


t=2  i=JV 


d()dry  E  r<r"V  E  fi{0)eiaeipeiy  =  0  (2-58) 
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The  square  bracket  in  the  equation  above  vanishes  in  accordance  to  Eq.  (2.57). 
If  we  take  into  account  Eqs.  (2.3)  and  (2.19)  and  consider  n  =  r2  =  r,  the 
second  order  momentum  equation  becomes 

dt2  Y  Pau'a  +  r(dtl)2  Y  pau'a 

<7=1  <7=1 

-  rkBT  Y  ( 2daV  •  u'  +  VVj  =  0  (2.59) 

<7=1 

When  adding  the  first  and  second  ordet  momentum  equations,  we  get 

<7=2  <7=2  <7=2 

dt  Y  P*U'a  +  T  (^)2  E  P°U'«  +  d«P  +  ^E  P<7<4 

<7=1  <7=1  <7=1 

-  rkBT  £  [29„(V  -u')  +  V2<]  =  0  (2.60) 

<7=1 

Thus,  the  momentum  equation  up  to  the  second  order  contains  the  first 
and  second  order  time  derivatives,  like  the  mass  equation  of  the  whole  fluid 
system  (2.45.  As  mentioned  previously,  the  presence  of  the  second  order  time 
derivative  is  a  characteristics  of  a  propagation  phenomenon  [18]. 

2.7  Pure  diffusion 

We  suppose  that  the  fluid  system  is  always  at  rest.  In  this  case,  the  equilib¬ 
rium  velocity  u'a  should  vanish  everywhere.  When 

<  =  0  (2.61) 

the  series  expansion  (2.21)  reduces  to 

/"<0)  =  ft™  =  wtn,  (2.62) 

Consequently,  the  first  order  mass  equations  (2.38)  becomes 

i=N  i=N 

0tl  Y  Win*  +  dP  S  Win*eil3  =  0  (2-63) 

£=0  i=0 


55 


In  accordance  to  Eqs.  (2.13)  and  (2.14),  we  have  the  following  equalities 


Y  wi  =  1 
1=0 
i=N 

Y  wieif3  =  0 

i=0 


(2.64) 


(2.65) 


which  are  valid  for  N  =  6,  as  well  as  for  N  =  8  (i.e.,  independent  of  the 
number  of  the  discrete  velocity  speeds  in  the  phase  space).  Thus,  the  first 
order  mass  equation  (2.63)  reduces  to 


dtl  na  =  0 


(2.66) 


Therefore,  from  the  first  order  Boltzmann  equation  (2.27)  we  get,  using  the 
expression  (2.62) 

fi{1)  =  -raWieaiadana  (2.67) 

which  may  be  introduced  in  the  second  order  mass  equation  (2.39)  to  give 

i=N  i=N  i=N 

dt2  Y  win*  -  dh  Y  TffWie?adana  -  dp  Y  Tcrwieiaebd<*no  =  0  (2-68) 

i=0  1=0  1=0 


Since 


i=N 

Y  Wieiaeb  =  Xa4Sa(3 

i=0 


the  second  order  equation  (2.68)  is  rewritten  as 

dt2na  =  raXaC2aVna 


(2.69) 


(2.70) 


After  multiplication  of  the  mass  equations  (2.63)  and  (2.70)  by  e  and  £2, 
respectively,  we  can  sum  them  together  to  recover  the  mass  equation  up  to 
second  order  with  respect  to  the  Knudsen  number 


dt  n,  =  T„x<7<£Vn0 


(2.71) 


which,  after  multiplication  with  ma,  becomes  identical  to  the  pure  diffusion 
equation 

dtPa  =  W2pa  (2.72) 
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where  the  diffusion  coefficient  is 


V  =  TaXaCl  =  Tt 


kBT 


m„ 


(2.73) 


If  we  remind  the  expression  of  the  viscosity  in  the  FDLB  model  [1] 


v  = 


kBT 


m„ 


(2.74) 


we  see  that  the  Schmidt  number  Sc  [12]  is  a  constant  equal  to  unity  in  the 
present  FDLB  model 

So  =  £  =  1  (2.75) 


2.8  Comparison  with  the  theory  of  Shan  and 
Doolen 

The  LGLB  model  previously  developed  by  Shan  and  Doolen  [13,  17]  was 
used  to  simulate  the  time  evolution  of  the  barycentric  velocity  profile  in  a 
diffusion  couple  [19].  Here  is  a  brief  outline  of  the  main  characteristics  of 
this  model: 

1.  All  components  have  the  same  equilibrium  velocity  u',  which  means 
that  the  hypothesis  expressed  by  Eq.  (2.31)  is  used. 

2.  The  equilibrium  velocity  u'  in  the  absence  of  external  forces  is  given 
by  Eq  (2.54). 

3.  The  equilibrium  distribution  functions  f*,eq  are  expressed  as  series  ex¬ 
pansions  in  the  equilibrium  velocity  u'  while  the  leading  order  distribu¬ 
tion  functions  are  expressed  as  series  expansions  in  the  local  fluid 
velocity  u,  (i.e.,  the  barycentric  velocity  -  see  [17],  as  well  as  [19]);  this 
means  that  Eq.  (2.29)  is  not  valid  in  the  model  of  Shan  and  Doolen; 
in  fact,  they  make  no  explicit  use  of  the  zero-th,  first  and  second  order 
Boltzmann  equations  (2.26  -  2.28)  as  separate  entities;  we  may  imagine 
that  the  starting  point  of  the  model  of  Shan  and  Doolen  is  the  sum  of 
the  zero-th  and  first  order  Boltzmann  equation 

9tJ-m  +  =  -  —  [/f(0)  +  f’m  -  (2.76) 

Ta 

which  allows  /f®  /  ff’eq,  and  so  on  [19]. 
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4.  The  series  expansion  for  /f’e9  and  which  is  valid  for  N  =  6, 
i.e.,  the  six  bit  model  (see  [19])  contains  the  parameter  da  which  may 
be  adjusted  to  allow  different  values  of  the  diffusivity;  this  “degree 
of  freedom”  does  not  exist  in  the  FDLB  model  since  the  very  recent 
approach  of  He  et  al.  [20,  21,  22]  clearly  establishes  the  value  da  =  0.5 
using  Gaussian  quadrature  formulae. 

5.  We  mention  here  that  Shan  and  Doolen  give  two  simulation  results 
(figures  1  and  2  in  their  paper  [13]);  these  figures  refer  to  equilibrium 
density  profiles  in  a  binary  mixture  which  seem  to  be  correct;  since  these 
results  refer  to  systems  at  equilibrium,  we  think  that  the  barycentric 
velocity  u,  as  well  as  the  equilibrium  velocity  u'  vanish  and  the  pure 
diffusion  case  is  recovered;  this  may  be  a  serious  argument  to  consider 
the  model  as  being  valid  only  for  stationary  (equilibrium)  cases,  when 
there  is  some  competition  between  diffusion  and  other  phenomena  gen¬ 
erated  by  external  forces,  while  the  model  itself  does  not  account  for 
the  true  dynamic  evolution  towards  the  equilibrium  (stationary)  case. 

6.  We  should  mention  also  the  fact  that  the  LGLB  model  of  Shan  and 
Doolen  inherits  the  main  disadvantage  of  the  LGLB  models,  which  is 
the  fact  that  the  particle  thermal  speeds  in  this  model  are  strictly  re¬ 
lated  to  the  lattice  spacing  (in  fact,  the  magnitude  of  these  speeds  is 
strictly  the  lattice  spacing  divided  by  the  time  step)  and  thus,  any 
LGLB  model  does  not  allow  different  thermal  speed  for  particles  car¬ 
rying  different  masses;  FDLB  models  allow  different  thermal  speeds 
when  the  particles  have  differents  masses  and  thus,  one  may  expect 
these  models  to  be  more  close  to  the  physical  reality. 

2.9  Simulation  results 

We  used  a  250  x  5  square  lattice  with  walls  placed  left  and  right.  The 
lattice  spacing  was  5x  =  Sy  =  10~4  m.  Periodic  conditions  were  used  at 
the  upper  and  bottom  boundaries.  The  left  half  of  the  lattice  was  initialized 
with  particles  belonging  to  species  a  =  1  while  the  right  half  was  initialized 
with  particles  belonging  to  species  a  =  2.  When  developing  the  computer 
code,  we  used  the  nine  bit  model  (N  =  8)  since  this  model  allowed  a  larger 
variety  of  finite  differences  schemes  to  be  tested  (see  Chapter  1). 
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Our  first  diffusion  code  used  the  Second  Order  Runge  Kutta  for  time  in¬ 
tegration  [1],  combined  with  an  Upwind  scheme  for  calculating  space  deriva¬ 
tives.  Later,  we  developed  a  code  based  on  the  First  Order  Upwind  Scheme 
(FU)  which  we  described  in  the  previous  chapter.  Both  the  Runge  Kutta 
code  (Appendix  B),  as  well  as  the  FU  code  (Appendix  A)  gave  similar  re¬ 
sults  which  are  described  below. 

The  initial  particle  number  density  was  set  to  unity  for  each  component. 
As  a  result,  the  initial  density  profiles  of  each  component,  as  well  as  the  total 
density  profile  were  always  similar  to  those  depicted  in  Figure  2.1,  where  we 
used  mi  =  1  amu  and  m2  =  0.9  amu.  1  amu  (atomic  mass  unit)  equals 
1.661  x  10"27  kg. 

The  present  simulations  were  done  with  the  same  value  of  the  relaxation 
time:  ra  =  r  =  10-8  s,  while  the  time  step  was  chosen  to  be  8t  =  10~9  s. 
The  thermal  velocities  of  each  species  of  particles  were  calculated  in  accor¬ 
dance  to  Eq.  (2.3)  with  T  =  300  K. 

Figure  2.2  shows  the  typical  time  evolution  of  the  total  density  profile. 
One  can  see  the  presence  of  two  kinks:  a  left  propagating  kink  and  a  right 
propagating  one.  Their  evolution  is  presented  in  Figures  2.3  and  2.4.  The 
presence  of  these  propagating  kinks  is  not  a  surprise  if  we  remind  the  mass 
equation  (2.45)  of  the  whole  fluid  system,  which  contains  the  second  order 
time  derivative.  As  mentioned  before,  the  presence  of  the  second  order  time 
derivative  is  a  characteristics  of  a  propagation  phenomena. 

Figure  2.5  shows  the  time  evolution  of  the  barycentric  velocity.  We  may 
see  the  two  kinks  which  form  and  propagate  laterally,  while  the  profile  ex¬ 
hibits  a  central  peak  which  remains  always  positive.  The  propagating  kinks 
are  reflected  by  the  lateral  walls  and  superpose  to  the  central  peak  at  some 
moments  during  their  propagation.  The  formation  of  these  propagating  kinks 
in  the  figures  showing  the  time  evolution  of  the  barycentric  velocity  may  be 
explained  by  the  presence  of  the  second  order  time  derivative  in  the  momen¬ 
tum  equation  (2.60). 

Figure  2.6  shows  the  time  evolution  of  the  mass  flux  of  component  2. 
This  flux  is  orientated  from  left  to  right,  as  expected,  and  does  not  change 
its  sign.  No  kinks  are  observed. 

Figure  2.7  shows  the  total  density  profile  when  the  two  species  of  par¬ 
ticles  have  the  same  mass.  This  profile  is  found  to  be  constant  during  the 
diffusion  process,  while  the  individual  profiles  of  components  1  and  2  evolve 
separately.  Figure  2.8  shows  the  corresponding  time  evolution  of  the  local 
density  of  component  1.  One  can  see  that  kinks  are  not  present  in  the  den- 
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Figure  2.1:  Initial  density  profiles  in  the  diffusion  couple:  a  -  component 
a  =  1;  b  -  component  a  =  2;  c  -  total  density  profile  p  —  p\  +  p2- 
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a:  1,000  time  steps 


b:  2,000  time  steps 


c:  50,000  time  steps 


Figure  2.2:  Time  evolution  of  the  total  density  profiles  in  the  diffusion  couple. 
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Figure  2.3:  Time  evolution  of  the  total  density  profiles  in  the  left  region  of 
the  diffusion  couple  (snapshots  are  taken  every  1000  time  steps). 
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Figure  2.3:  (continued)  Time  evolution  of  the  total  density  profiles  in  the  left 
region  of  the  diffusion  couple  (snapshots  are  taken  every  1000  time  steps). 


Figure  2.4:  Time  evolution  of  the  total  density  profiles  in  the  right  region  of 
the  diffusion  couple  (snapshots  are  taken  every  1000  time  steps). 
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Figure  2.4:  (continued)  Time  evolution  of  the  total  density  profiles  in  the 
right  region  of  the  diffusion  couple  (snapshots  are  taken  every  1000  time 
steps). 


sity  profiles  when  masses  are  equals,  which  means  that  we  are  dealing  with 
a  pure  diffusion  process. 

To  compare  the  FDLB  simulations  wit  the  former  LGLB  simulations  [19], 
we  slightly  modified  the  corresponding  computer  code  and  set  the  same  ther¬ 
mal  velocity  for  each  species  of  particles,  even  if  the  masses  are  different. 
For  convenience,  we  used  mi  =  1  amu,  m2  =  0.9  amu,  while  C\  was  given 
by  Eq.  (2.3)  and  c2  was  forced  to  be  equal  to  c\.  Figure  2.9  shows  the  time 
evolution  of  the  barycentric  velocity  in  the  diffusion  couple,  which  is  very 
similar  to  the  LGLB  results  reported  in  [19].  In  this  case,  the  barycentric 
velocity  changes  its  sign  alternatively,  due  to  strong  oscillations  which  are 
present  in  the  system. 

The  unphysical  behavior  of  the  barycentric  velocity  reported  in  figure  2.9, 
as  well  as  the  simulation  results  reported  in  [19]  provide  a  very  strong  reason 
for  a  careful  handling  of  the  former  LGLB  models  when  dealing  with  mul¬ 
ticomponent  (e.g.,  binary)  systems  whose  particles  carry  different  masses. 
All  the  present  Lattice  Boltzmann  litterature  dealing  with  multicomponent 
systems,  e.g.,  the  approaches  in  [11,  14,  18,  23,  24,  25,  26,  27,  28,  29,  30],  ig¬ 
nores  completely  the  connection  between  the  lattice  spacing  and  the  thermal 
speeds  of  particles  and  do  not  report  unphysical  effects  of  the  Lattice  Gas  like 
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Lattice  Boltzmann  scheme.  The  unphysical  oscillations  which  are  observed 
when  trying  to  simulate  the  behavior  of  a  diffusion  couple  in  microgravity  en¬ 
vironment  with  the  LGLB  model  revealed  a  severe  limitation  of  these  models 
to  those  multicomponent  systems  whose  particles  have  identical  mass  (as  well 
as  the  same  thermal  velocity).  Finite  Difference  Lattice  Boltzmann  Models 
(FDLB)  or  Interpolation  Supplemented  Lattice  Boltzmann  Models  (ISLB) 
should  be  used  in  order  to  allow  different  thermal  velocities  of  particles  in 
the  system  (i.e.,  different  values  of  the  Courant  -  Friedrichs  -  Lewy  number 
CFL).  In  this  respect,  the  absence  of  the  oscillations  of  the  barycentric  veloc¬ 
ity  in  the  diffusion  couple,  as  reported  in  Figure  2.5,  are  very  encouraging. 
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t  -  2000 


Figure  2.5:  Time  evolution  of  the  barycentric  velocity  profile  (snapshots  are 
taken  every  1000  time  steps). 
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Figure  2.5:  (continued)  Time  evolution  of  the  barycentric  velocity  profile 
(snapshots  are  taken  every  1000  time  steps). 
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Figure  2.6:  Time  evolution  of  mass  flux  of  component  2. 
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Figure  2.7:  Total  density  profiles  in  the  diffusion  couple  when  the  two  species 
of  particles  have  equall  mass:  this  profile  remains  constant  during  the  diffu¬ 
sion  process. 


Figure  2.8:  Time  evolution  of  the  density  profiles  of  component  1  when  the 
two  species  of  particles  have  equall  mass:  no  kinks  are  observed. 
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Figure  2.9:  Time  evolution  of  the  barycentric  velocity  profile  (snapshots  are 
taken  every  5000  time  steps,  except  the  second  one,  which  is  taken  after  1000 
time  steps)  when  the  two  thermal  velocities  are  set  equal  (ci  =  c2). 
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Chapter  3 
Sessile  drop 


3.1  Drop  shape 

A  computer  code  based  on  a  lattice  Boltzmann  model  for  liquid  -  vapor  sys¬ 
tems  was  already  described  in  the  preliminary  report  [31].  This  code  was 
used  to  investigate  the  equilibrium  shape  of  a  liquid  drop  on  a  horizontal 
surface  and  some  preliminary  results  were  presented  in  [31].  These  prelimi¬ 
nary  simulations  were  done  on  a  256  x  64  lattice  because  of  limitations  in 
the  computer  resources  available  at  that  time.  In  the  present  report  we  will 
refer  to  subsequent  simulations  made  on  larger  lattices,  which  were  done  on  a 
Pentium  II  processor  with  256  MB  RAM  working  under  the  UNIX  FreeBSD 
operating  system.  The  code  is  given  in  Appendix  C. 

Following  constant  values  of  the  main  simulation  parameters  were  used, 
as  in  [31] 

•  mean  system  density  pm  =  3.5 

•  conventional  temperature  T  =  0.550 

•  relaxation  time  r  =  0.8 

To  investigate  the  shape  of  a  sessile  drop  on  a  horizontal  surface  subjected 
to  a  gravitational  field,  we  used  a  constant  value  of  the  parameter  k,  =  0.01 
of  the  LB  model,  which  defines  the  surface  tension,  as  discussed  in  [31].  We 
considered  the  two  different  cases,  when  the  liquid  is  wetting  the  horizontal 
surface  or  not.  For  the  first  case,  we  adopted  the  value 

A-  =  0.01  (3.1) 

uXji 
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of  the  normal  derivative  of  the  chemical  potential  on  the  lower  wall  surface. 
This  value  of  the  derivative  accounts  for  the  wetting  interaction  between  the 
liquid  drop  and  the  wall.  The  second  case  is  characterized  by 

4r~  =  -o-oi  <3-2) 

OXfi 

which  accounts  for  a  non  -  wetting  wall. 

Figure  3.1  shows  the  time  evolution  of  a  sessile  drop  which  wets  the 
bottom  wall  and  is  subjected  to  an  acceleration  directed  downwards  ( g  = 
-0.00001).  The  drop  spreads  across  the  wall  surface.  Because  the  lattice 
was  not  large  enough,  the  drop  ends  join  laterally  because  of  the  periodic 
boundary  conditions  we  used  in  the  X  direction,  so  that  a  flat  liquid  surface 
is  recovered  in  the  final  state. 

Figure  3.2  shows  the  equilibrium  shape  of  a  nonwetting  drop  subjected 
to  different  values  of  the  gravitational  acceleration.  One  can  see  that  the 
drop  becomes  more  and  more  flattened  when  the  value  of  the  gravitational 
acceleration  increases. 


3.2  Contact  angle 

We  used  a  1024  x  192  lattice.  The  necessary  CPU  time  to  perform  100,000 
time  steps  was  approx.  72  hours  for  each  data  set.  Several  values  of  the 
surface  tension  coefficient  k  were  considered,  while  the  normal  derivative  of 
the  chemical  potential  on  the  lower  wall  surface  was  =  0.01.  To  account 
only  for  the  effect  of  surface  tension  on  the  contact  angle,  no  gravitational 
acceleration  was  considered  (g  =  0). 

Figures  3.3  -  3.6  show  the  time  evolution  of  the  lattice  state  for  different 
values  of  the  surface  tension  coefficient  k.  When  the  value  of  the  surface 
coefficient  is  large  enough  (k  =  0.01),  the  drop  spreads  along  the  whole 
lattice  domain,  so  that  the  contact  angle  cannot  be  determined  because  of 
the  periodic  boundary  conditions  on  the  left  and  right  margin  of  the  lattice 
(Figure  3.3).  For  smaller  values  of  the  surface  tension,  coefficient,  when  the 
drop  does  not  spread  along  the  whole  bottom  wall,  the  contact  angle  a  may 
be  determined  using  the  following  procedure  (Figure  3.7),  which  takes  into 
account  the  fact  that  the  upper  drop  border  becomes  a  circular  arc  when  the 
equilibrium  configuration  is  reached  (i.,e.,  when  the  number  of  time  steps  is 
large  enough).  The  drop  height  is  denoted  h,  while  its  width  is  denoted  w. 
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t  =  2,000  time  steps 


Figure  3.1:  Time  evolution  of  the  shape  of  a  wetting  drop. 
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t  =  3, 000  time  steps 


t  =  4, 000  time  steps 


t  —  5, 000  time  steps 


Figure  3.1:  (continued)  Time  evolution  of  the  shape  of  a  wetting  drop. 
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t  =  10, 000  time  steps 


Figure  3.1:  ( continued)  Time  evolution  of  the  shape  of  a  wetting  drop. 
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Figure  3.1:  (continued)  Time  evolution  of  the  shape  of  a  wetting  drop. 


If  R  is  the  arc  radius,  we  have 


h 

=  i?(l  —  cos  cc) 

(3.3) 

w 

=  2 R  sin  a 

(3.4) 

h 

1  —  cos  a 

1  a 

(3.5) 

—  tan  — 

w 

2  sin  a 

2  2 

cos  a 

(  2  h 

—  2  arctan  - 

) 

(3.6) 

\  w 

Table  3.1  shows  the  values  of  the  contact  angles  for  the  cases  shown  in 
figures  3.4  -  3.6. 
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g  =  -0.000001 


g  =  -0.000020 


g  =  -0.000050 


Figure  3.2:  Shape  of  a  non  -  wetting  drop  subjected  to  different  values  of  the 
gravitational  field. 
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initial  state  ( t  =  0) 


t  =  1, 000  time  steps 


t  =  5, 000  time  steps 


Figure  3.3:  Drop  shape  during  simulation  with  k,  =  0.0100. 
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t  =  20, 000  time  steps 


t  =  30, 000  time  steps 


t  =  32,000  time  steps 


t  =  33, 000  time  steps 


Figure  3.3:  (continued)  Drop  shape  during  simulation  with  k  =  0.0100. 


82 


initial  state  ( t  =  0) 


t  —  1,000  time  steps 


t  =  5, 000  time  steps 


t  =  10, 000  time  steps 


Figure  3.4:  Drop  shape  during  simulation  with  k  =  0.0094. 


initial  state  ( t  =  0) 


t  =  1,000  time  steps 


t  =  5, 000  time  steps 


Figure  3.5:  Drop  shape  during  simulation  with  k  =  0.0090. 
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t  =  100, 000  time  steps 


Figure  3.5:  (continued)  Drop  shape  during  simulation  with  k  =  0.0090. 
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Figure  3.6:  (continued)  Drop  shape  during  simulation  with  k  =  0.0088. 
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Table  3.1:  Influence  of  the  surface  tension  coefficient  k  on  the  contact  angle. 


K 

cos  a 

0.0094 

0.9097 

0.0090 

0.7625 

0.0088 

0.6671 

Figure  3.7:  Determination  of  the  value  of  the  contact  angle  a. 
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